% This file performs a Monte Carlo exercise for different values of r and
% lambda that are being assigned in the main file
% 'FIGs_D2_D3_Monte_Carlo_Varying_r_lambda.m', creating Figure D2 and D3
% presented in Appendix D of the paper. 

%% Generate the True DGP. 

A           =   [0.7 0.1;
                0.2 0.8];
B           =   [1 0;0.2 1];

H           =   9;                      % Number of horizons.
IRF         =   zeros(2,H);
for jx = 1:H
    IRF(:,jx)   =   A^(jx-1)*B(:,1);
end

%% Generate Fake Data.  

N           =   1000;                   % Number of Monte Carlo Simulations.
TT          =   200;                    % Number of periods.
burn        =   100;                    % Number of periods to burn-in. 
T           =   TT + burn;              % Total number of periods to generate.
p           =   1;                      % Number of lags.

randn('seed',727);
e           =   1*randn(2,T,N);
Zsim        =   zeros(2,T,N);
IRFmc       =   zeros(2,H,N);

for jx = 1:N
    Zsim(:,1,jx)    =   B*e(:,1,jx);
    for ix = 2:T
        Zsim(:,ix,jx)    =   A*Zsim(:,ix-1,jx) + B*e(:,ix,jx);
    end
    [AA,SIGMA,U,V]  =   olsvarc(Zsim(:,burn+1:end,jx)',p);
    B2              =   chol(SIGMA(1:2,1:2));
    B2              =   B2';
    [IRFM]          =   IRFVAR(AA,B2,p,H-1);
    IRFmc(:,:,jx)   =   IRFM(1:2,:);
end

% Compute average across runs:

IRFmca      =   zeros(2,H);

for kx = 1:2
    for ix = 1:H
        IRFmca(kx,ix)   =   mean(IRFmc(kx,ix,:));
    end
end

irflp       =   zeros(N,H);
irfslp      =   zeros(N,H);

% Now do a local projection of Y on X:
for jx = 1:N
    x           =   Zsim(1,burn+1+p:end,jx)';
    y           =   Zsim(2,burn+1+p:end,jx)';
    w           =   [Zsim(1,burn+1:end-1,jx)' Zsim(2,burn+1:end-1,jx)'];
    h1          =   1;
    lp          =   locproj_var(y,x,w,h1,H,'reg'); 
    irflp(jx,:) =   lp.IR(2:H+1)'; 
    slp         =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
    % Extract the IRF from local projection: 
    irfslp(jx,:) =  slp.IR(2:H+1)';
end

irflpm      =   zeros(1,H);
irfslpm     =   zeros(1,H);
for jx = 1:H
    irflpm(1,jx) =  mean(irflp(:,jx));
    irfslpm(1,jx)=  mean(irfslp(:,jx));
end